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ABSTRACT 



A detector system locates a feature in an input signal. The 
input signal is applied to a multiscale filler bank which 
produces a corresponding set of filtered signals that are 
applied to a multiscale detector. A scale space tree table is 
produced by the detector which indicates the location in 
each filtered signal where the feature is possibly located. A 
MAP correlator examines the scale space tree table to 
determine which indications therein are feature locations. 
One-dimensional and two-dimensional embodiments are 
described. In addition, the design framework is disclosed for 
extending a one scale filter to a multiscale filter, for design- 
ing a filter bank from the multiscale filter and for designing 
the MAP correlator dependent on the feature being located 
and the chosen filter bank. 

6 Claims, 7 Drawing Sheets 
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MULTISCALE FEATURE DETECTOR USING in the image. Noise is a very significant factor, for example, 

FILTER BANKS with medical images and the present invention is particularly 

useful in such applications. 

CROSS-REFERENCE TO RELATED Another object of the present invention is to match the 
APPLICATION 5 MAP correlator with the mnltiscale filter so that optimal 

This application is based on ProvisioDal Application No. performance is achieved. The model response and the preset 

60/027,422 filed on Jul. 29 1996 threshold used by the MAP correlator to determine the 

"Tliis invention was made with United States Govermnent P^.'^"^ ^ feature are precisely determined by the oper- 

support awarded by the NaUonal Science Foundation: NSF ^^""^ parameters of the chosen multiscale filter bank. 

Grant No.: MIP-9501589. The United States Government Th^ foregoing and other objects and advantages of the 

has certain rights in this invention." invention will appear from the following description. In the 

description, reference is made to the accompanying draw- 

BACKGROUND OF THE INVENTION ings which fonn a part hereof, and in which there is shown 

™_ cij • ...... rr . by way of illustration a preferred embodiment of the in ven- 

Ihe field of the invenUon is the dctecUon of features m 15 ^^^^ embodiment does not necessarily represent the 

images, and particularly the detection of features such as invention, however, and reference is made 

edges m two-dimensional images. ^j^^^^^^^ ^^^.^^ ^^^^-^ interpreting the scope of 

There are numerous applications in which images are the invention, 
analyzed to locate specific features. For example, an object 

may be scanned with a camera and the edges of the object BRIEF DESCRIPTION OF THE DRAWINGS 

may be detected and used to locate the object and determine ^ • j- 

its orientation. Edge detection is also used in estabhshing the ^ ^ ^mgr^m of the elements of a feature detection 

boundaries of an object so that it may be identified using system, 

pattern recognition methods, or measured. For example, an 2 is a graphic representation of discrete multiscale 

ultrasound image of the human heart may be analyzed by filters at different sampling rates and their corresponding 

first locating the edges of the various heart chambers and frequency transforms derived from the Haar wavelet; 

then calculating4he sizes of those chambers and the thick- FIG. 3 is a diagram of a multiscale filter; 

nesses of the chamber walls. FIG. 4 is a schematic representation of the relationship 

There are many approaches used to detect features in between a multiscale filter and a corresponding multiscale 

images. One approach is to fit a model of the feature to the filter bank; 

image. The model is moved around the image and where FIG. 5 is a graphic representation of a sampling lattice 

there is a good match the feature location is indicated. Other produced by a decimator D used in a 2D filter bank of the 

methods use probabilistic measures, and most methods use present invention; 

filtering tcdnuques which attempt to remove noise without g ^ ^.i^aiic^i block diagram of the preferred 

obscurmg the feature being detected Filter operators such as embodiment of the invented feature detection system; 
the Roberts operator, Haar operator, Sobel operator. Canny s 

detector and wavelet transform are used to enhance images ^ ^ schematic representation of an mdicator 

as a first step in feature detection, but the effectiveness of response table which is a data structure employed in the 

such filtering is 'determined by the choice of the filter scale. ^^^^^"^ of FIG. 6; 

The choice of the filter scale may be arbitrary, or it may be 8 is a schematic representation of a scale space tree 

based on results from prior applications of a similar nature. ^^^^^ which is a data structure employed in the system of 

If the filter scale is large, filtering is intense and noise is ^'y 

smoothedoutatthecxpcnseof lost localization of a detected FIG. 9 is a flow chart of the functions performed by a 

feature or even failure to detect a feature. If the scale is feature generator which forms part of the system of FIG. 6; 

small, on the other hand, less filtering is used and the results and 

may be corrupted by noise which produces false alarms. fig. 10 is a flow chart of the functions performed by the 

SUMMARY OF THE INVENTION correlation which forms part of the system of FIG. 6. 

Hie present invention is a method and apparatus for 50 GENERAL DESCRIPTION OF THE INVENTION 

detecting features in an acquired image, and particularly, it The present invention is a system for detecting and 

includes: a multiscale filter bank which applies a filter locating features in an input signal. Referring to HG. 1 it 

operator to the image data to produce a plurality of filtered includes a multiscale filter 10 which receives the input signal 

outputs at a corresponding plurality of scales; a detector for 12 and filters the signal to produce a set of filtered output 
each filtered output which indicates the locations of a 55 signals 14a-c. The filtered output signals 14fl-c represent 

selected feature therein, and a MAP con-elator which cor- the same input signal 12, but they are filtered at different 

relates the filtered output values at each indicated feature "scales". ITie lower scale signal 14a is the highest resolution 

location with a model response and produces an output signal and it has undergone the least filtering. The highest 

signal indicative of the presence of the feature at the location scale signal 14c has been filtered the most and it is the lowest 

if the correlation exceeds a preset threshold. resolution signal. While ID signals are shown in FIG. 1, 2D 

A general object of the invention is to provide a feature signals such as images and 3D signals such as image 

detection system which is robust. Nearly all state-of-the art sequences (video) may also be processed, 

feature detectors work satisfactorily when the image is not The multiscale filtered signals 14a-c are apphed to a 

corrupted with noise. multiscale detector 16. The detector 16 provides an indica- 

However, using the multiscale filter bank and the map 65 tion of the location in the filtered signals 14fl-c of the 

correlator of the present invention enables the same features desired feature. "Features" to be detected may include, for 

to be accurately located even when there is substantial noise example, step edges, roof edges, and pulse edges, and these 
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may be "indicated" by maxima, zero-crossiag, or a constant 

interval in the filtered signals. For example, maxima and /--i; -i<r<o 

zero-crossings are two common indicators for detecting ^(0 = J l; 0 </< i 

edges and the choice of indicator will depend on the feature I q otherwise, 
to be located and the type of filter used. For example, when S 
the feature is a step edge, a maxima detector may be paired 

with a filter such as that described by Canny, "first diflference Samples of this wavelet give a type 1 filter that is 

operator, first derivative of Gaussian" or a quadratic spline symmetric with even length. Assume that the largest support 

filter. On the other hand, a zero-crossing detector may be S^^^r has length 16. The Haar wavelet is uniformly sampled 

paired with a Laplacian of Gaussian (LOG) filter. Regardless lO samples as shown in FIG. 2 at rate 4 (rate and 

of the choice, the multiscale detector 16 produces samples of ^^^^ mterchangeably here). The smaller support 

the filtered signals I4a-c at each indicated location at ^^^.^^"^^^ ^'^^''.'^1^'°!'^'°^ f • '"'^1.^^ nZ^l 

corresponding outputs IHa-c, Tlie multiscale detector 16 in 'T^^\ J' i f t f .^^^^1^*^°° ^oif) 

essence idenUfies locations in the filtered signals 14a-c ^^^^^ ^^^^^ ^^^^^^ of samples, one needs to 
. r * u * J 1 rTu • , examme the frequency responses. We observe the following 

where the feature may be present and masks out the remain- is ~ . .. . r« ^ , % 

der of the signal. TUis step significantly reduces the pro- f""^^ f '^f'f"' ^'^P'""^ J^, ^"^^ S""' °^ 

, t w 1 IX AH 1 . iA u- u length two at rate 1, passes components at frequency k. 

cessing burden on a multiscale MAP conelator 20 which tt j cfl c^ r . . ^ j . 

tk^ ^^t^™;„of;«« if tu^ V ««o^«f However, the second filter of length four at rate 2 does not. 
makes the final celermination if the feature is present. r™. • ■ A . • r • • . l ^ . ^zi. ■ 
^ , . , , ^« ^ , ^ , , . This is a first sign of inconsistency where the first filter is not 
The muluscale MAP correlator 20 makes the final deci- bandlimited. In addition, the z-transforms of both filters do 
sion of where the feature occurs in the signal 12 and 20 seem to correspond to filters related by sampling. This is 
produces an output 22. As will be described in more detail another critical sign of inconsistency. On the other hand, 
below, It makes this determinaUon by correlating the detcc- comparing the second filter to the third, and so on, consis- 
tor output signals ISa-^ with a reference vector that repre- jency is observed in the responses. The conclusion is that the 
sents the expected response when the feature is present. This gj^^Hest filter, where consistency starts to show, must be of 
reference vector as well as the threshold value of the 25 ^^^^^^ ^^^^ ^^^^ Detection of features using a multi- 
correlation required to indicate the presence of the feature is 3^.^]^ ^jter is reinforced from the consistency of response at 
precisely defined. In a ID system the output 22 indicates the ^^^^^^^ ^ ^^e filter does not provide this property, then 
times at which the feature occurs in the signal 12, whereas ^ ^^i^gble for multiscale feature detection according to 
in a 2D system the output 22 is a bit map showing where in j^e present invention. Other filters that may be employed in 
the 20 mput signal 12 the feature is located. 30 ^^^^^^ invention include: 

The construction of the multiscale filter for this system Laplacian of Gaussian (LOG), where the smallest filter 

must meet a number of design criteria. The highest rcsolu- has length 17; 

tion filter is the smallest filter in the multiscale set, and it is fi^st derivative of Gaussian (FOG), where the smallest 

the most critical for localization of the feature. To determine gj^^j. ^ length of 9* 

the highest resolution scale, one needs to "examine the triangular waveform, where the smallest Mter has a length 

Fourier transforms ot the desired set of filters. The filters in 5. 

the multiscale set have to preserve their shape and ^ ' i_ r 1 .1 n . 01. 1 1 . 

j*u- .-.u- sawtooth waveform, where the smallest filter has a length 

symmetry, and this should be apparent in their responses. of 9 

This is a necessary condition for the extension of one -scale ™ ^i' * - ^ ^ . .1 ^i. t . 

, !*• 1 o* *• f *u 1 * u f 40 The niters in this set of acceptable filters are denoted 

filters to multiscale. Startmg from the least number of . ..^/\^ ^ ■ 1 . ^ - 

1 , f 1 j *u f j j 1 4u herein by G J z). These filters can be implemented id parallel 

samples, two tor even-length niters and three lor odd-length , • f-ry- 'j l . a/ - . . . • . 

^ J ^ P CI. • * * J as shown in FIG. 3, but a more emcient structure is to 

filters, the z transforms of the filters m the set are computed. . , , tiu ■ ci. u 1 u • r^r^ a 

™ , ' 1 . . V implement these filters in a filter bank as shown in FIG. 4. 

The lowest scale is then chosen where: ,. c / \ 

^ ^ ^ ... For some choices of G^z) an exact filter bank iraplemen- 

1. The frequency transform starts showing consistency at ^^^^^^ ^ possible, and in other cases an approximation is 
the higher scales. This is verified by checking the ^^^^ ^1 5^ discussed below. 

sampling relation between successive scales. The -j^^ concepts discussed above may be extended to cover 

"pulse" waveform is an example of a one-scale filter 2D domain. It is assumed that a one-scale flUer has been 

where consistency is not found, and accordingly the j^rived for optimal detection of the desired feature. The 

samples are unsuitable for multiscale feature detection. ^y^^^ either separable or non-separable. For example, 

2. The filter's response must be bandlimited. This is separable and non-separable Roberts and Sobel operators 
necessary since the multiscale filters will be imple- may be used to detect vertical and diagonal edges respec- 
mented using analysis filter bank. The cascade combi- lively. The size of such filter operators is chosen to cover the 
nation of filters becomes equivalent to reconstructing a largest scale of interest. 

larger support filter from a smaller filter. This is a 55 This 2D filter is then used to generate a set of multiscale 

special case of reconstruction, and aliasing must be 2D filters by applying different sampling rates. The samp hng 

avoided. lattice can be characterized by a 2x2 nonsingular integer 

3. The highest resolution filter corresponds to the filter matrix. In ID, the relation between a signal x (n) and its 
with the smallest support. Therefore, this filter's decimated version x/n) by a factor M is given by: x,(n)«x 
response satisfies the above two constraints, and has its 60 (Mn). In 2D, the decimation factor is given by the determi- 
band width closest to jt. This results in a maximally nant of the decimation matrix D. The relation between the 
decimated filter, and higher scales are oversamples. 2D signal and its decimated version is given by: x^(n)=x 

The details of this procedure are now described for the (Dn). 

Haar filter, and similar methods are applied to other filter Two conditions are imposed on the decimation matrix D: 
operators without detaDing the steps. 65 In the ID approach described above, the decimation rate 

The sampUng rules described above are applied to the is 2. This choice is common in multiresolution analysis 

Haar wavelet as follows: for dyadic scales, and can be easily implemented using 
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two channel filter banks. By analogy to the ID case, the An exact solution for Ho(z) exists and is given by: 
2D rate should achieve a sampling density as power of 

four. This results in the requirement that the determi- H ii)- ^<^o(z) 0) 

nant of D is 4. 
However, the choice is not unique. This leads to the 

second condition which limits the alternatives. 
In the development of the detector, it is assumed that an rooXs of 

isolated feature indicator persists at all the scales. This 

implies that the decimator can not alter the shape of the Ci (zA + Ci f-zM 

filter. ^ ' ^ ^ 

Combining these two constraints leads to the unique deci- 
mation matrix given by: are also roots of Go(z). 

This is demonstrated as follows. Since Gq can be obtained 
[2 01 - ^5 by subsampling Gj, then 



10 2 



Co(z)=i[c,(zUG,(-z5)]. 



The generated sampling lattice resulting from this decima- 
tion matrix is shown in FIG. 5. Note that this lattice does not addition G Cz) H fz^^lH fz) Combinin the two 
necessarily lead to separable filters. For example, 2D filters , .' . , t*° (n\ 

with symmetry along diagonals can be designed to detect , , ./-^r 

r * ^ AC ji 1 TT.' • 4 i-i *L rhc above result may not always be satisfied. In addition, 

features at 45 degree angles. This is not possible with the tt / \ • * * j . u t l rr • 

X ^ .r^T^/-. .-. . Ho(z) IS not guaranteed to be linear phase. The Haar filter is 

separable processing approach of ID filters which restncts ^ ,^ ^^^^^ requirement is satisfied, and a linear 
the detection to honzontal and vertical directions only. ^^^^^ ^ ^^^^ 

Similarly, other non-separable filters can also be derived 

using the 2D method. ^^^^^ ^ _j ^^-z ^^-3 

With this choice of sampling lattice, a one-scale 2D filter 
generates a set of 2D discrete filters, denoted by Gj- for each 30 OiU) = -l-z~ -z~ -z' +z +z~ +z~ -t-z' 

scale i. These filters are related by the decimation matrix D: ; ! • 

G,^z)=Gy^i(z^. This decimation relation lends itself to a 
natural implementation as a filter bank. To derive this set, the 

one-scale 2D filter provides support large enough to cover and the exact solution is: 
the largest scale of observation and the rest of the filters are ^5 

successively derived by decimation on the lattice of FIG. 5. Ho(r)-i+2/*+r-2 and H,(r)^i-z-Vr-2+r-^ (3) 

As with the ID filter, the multiscale 2D filter G, may be ^ . , r .1.1. 

, . 1 . 11 1 i_ . i_ 1 ^- • \ The Haar filter is one example of an exact solution, but 

implemented m parallel, but a better solution is a tree . -ui r 1 •* ■ 1* * * j 

^ ^, , \ , , , t t others are possible. In general, it is difficult to extend 

structure filter bank as will be described below. i ci* * ™ u- i tin f • *l u i 

40 one-scale filters to multiscale filters satisfying the filter bank 

In designing the ID filter bank, the relation between the constraints. However, the reverse process is straightforward, 

parallel and serial muhiscale filter implementation in FIG. 4 A one-scale filter can be easily drawn from a filter bank by 

is discussed first. Given the desired set of filters G^z), one a simple collapsing of the cascaded filters. This is achieved 

needs to find the FIR filters Ho(z) and HJ^z) such that: by convolution of frequency multiplication. 

45 In cases, where the requirement in equation (2) is not 

Oo(z) = //i(z) (1) satisfied, an approximate solution based on least-squares 

2 approximation can be made. The desired filter Ho(z) should 
Oiiz) = Hoiz )Hi{z) chosen to minimize the error function: 



C?2(z) = //o(i^)Wo(2^)Wl(z) 



yriflo) = Y,^i/U)\\gi-g,.i*hoi\\ 



Hi(z) must match Go(z), and since application restricts the 

choice of Go(z), Hi(z) must follow. Tliis gives the first part where ho, is the upsampled version of ho, r is the number of 
of the filter bank, and to provide symmetry in the response, scales used in the analysis, and L,- is the length of g,. 
it needs to be a Jinear phase FIR. For the other part Ho(z), The following can. also be proven. Let N+1 be the size of 
several restrictions can be drawn by examining the structure. filter hp, E be an extension matrix of size (2N+l)x(N+l), U, 
It should be noted that a cascade of filters, as shown in FIG. be an upsampling matrix of size (L,-^j-L,.-l)x(L;^i-L,-)/2, 
4, is equivalent to one of the G^z) filters. The properties of and G,- the convolution matrix of size L,.^iXL.-. They are 
the equivalent filter at any scale should not be altered by the 
cascade leading to the next scale. For example, if Hi(z) is 
even-length and Ho(z^) is even-length, then the combination 
would not have an even-length and this would contradict the 
desired Gj(z). Therefore, Ho(z) needs to be symmetric of 65 
odd length to preserve characteristics, such as symmetry and 
anti-symmetry, across scales. 
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Furthermore, 



(4) 



V^hg) is minimum if ho=BoC"^a'. 

The above procedure may be used to design a filler bank 
for ID feature detectioQ based on the set of multiscale filters 
that are derived from a one -scale filter. This filter bank can 
also be used for 2D applications by applying separable 
processing, however it is not suitable for 2D applications 
where non-separable filters are required. In addition, 2D 
filters have different support regions than ID filters, and 
therefore the responses may be different. 

Conceptually, the formulation of a 2D filter bank is 
similar to the one dimensional case discussed above. Deci- 
mation matrices replace integer decimation, and matrices 
replace vectors in the signal representation. The set of filters 
G,-(z) is generated from a given one-scale 2D filter as 
explained in the previous section. The index i corresponds to 
the scale of observation, and it reflects the support size of the 
filter G-. One n^eds to design Ho(z) and H_i(z) such that: 



CoU) = //iU) 



IS 



(5) 



8 



where the hypothesis is satisfied, and an exact FIR filter 
solution can be found. 

In cases where the hypothesis fails, an alternative solution 
is proposed based on least-square approximation. The 
desired filter Ho(to) needs to minimize the error function: 



VriHo) 



\\\Gi-Gi,rHoir 



where Hq,- is the upsampled version of Hq, r is the number 
of scales used in the analysis, L,.=M,N,- is the size of the filter 
0;, and ** denotes 2D convolution. 

It can be shown that V/Hq) is minimum if Hq of size 
MxN is given by; 



Wo 



[Qvj^ where ho - EC a 



(6b)- 



20 



with 



25 



and 



Ol N 



30 



Moreover, the phase and length of the desired filters need to 
satisfy certain constraints as follows. The length of is the 
same as that of Gq, and the length of Kq can be derived from 
the lengths of Cq and Gj. Since symmetry and shape are 
preserved across the fillers G^, the same is expected when 
upsampled versions of Hq are cascaded with H^. This 
accordingly constrains Hq to be symmetric. 
An exact solution Ho(a)) exists and is given by 



//o(£)'(j) = 



l-^, Z GaD'^o>-'hik))\ 



C6a) 



if the roots of the denominator are also roots of the numera- 
tor. J(D) represents the sampling density of the decimeter D, 
and N(D^ is the set of integer vectors that lie inside the 
fundamental parallelepiped of D^ This can be demonstrated 
as follows. Gq is the resulting filter when Gj is cascaded with 
the decimator D. In other words, it is the decimated filter 
version of G, and 



On the other hand, from the equivalence of the filterbank 
structure to the parallel implementation Gj(ci))=Ho(D^tD)Hi 
(co). Combining these two equations yields equation (5). 

The above hypothesis may not always be satisfied and 
hence the exact solution of equation (5) may not work. In 
addition, Ho(a)) is not guaranteed to satisfy the desired 
symmetry constraints. The 2D Haar filter is an example 
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Qi is the 2D convolution matrix corresponding to the 2D 
fiher Gf, U,- is the upsampling matrix corresponding to the 
rectangular lattice decimator, 



45 



E is the matrix used to impose the symmetry constraints on 
the filter, where ho=»Ehoo with hoo=C~^ a. 

The dimensions of these matrices are as foUows: v„ is 

- g, is 



50 Nxl, N„ is MNxl hpQ is MqqXNqQ, 

M,N;Xl, a, and a^ a are MooNqoxI, 



is(M,,j+l-M,)(N,,,+l-N>(M,+l-M,_,)(N,+ 
1-N,.,), and B, is (M,,,+l-M,)(N,,i+l-N,)x(MooNoo). 
In summary, a 2D filter bank is designed by the following 
55 steps. 

Choose or design a 2D filter with a desired indicator 

response for the feature of interest. 
Generate a set of 2D filters sampled on a rectangular 

lattice as described above. 
Apply equations (5) and (6a) or (6i>) described above to 
determine the equivalent filter bank. 
These steps were followed to implement a 2D filter using a 
2D vertical Haar filter bank, using a horizontal Haar filter 
65 bank, a 2D vertical Sobel filter bank, and a 2D diagonal 
Sobel filter bank. The results of this implementation are 
summarized in Table 1. 
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TABLE 1 



Vertical Haar 


Ho 


row 


1 


1 2 1 






row 


2 


0 0 0 






row 


3 


1 2 1 




Hi 


row 


1 


-1-111 






row 


2 


-1-111 






row 


3 


-1-111 






row 


4 


-1-111 


Horizontal Haar 


Ho 


row 


1 


1 0 1 






row 


2 


20 2 






row 


3 


1 0 1 






row 


1 


-1 -1 -1 -1 






row 


2 


-1 -1 -1 -1 






row 


3 


1111 






row 


4 


1111 


Vertical Sob el 


Ho 


row 


1 


0.6956 1.3S91 0.956 






row 


2 


0.0242 0.0362 0.0242 






row 


3 


0.6956 1.3591 0.6956 




Hi 


row 


1 


-1-10 11 






row 


2 


-1-10 11 


• 




row 


3 


-2 -2 0 2 2 






row 


4 


-1-10 11 






row 


5 


-1-10 11 


Horizontal Sobel 


Ho 


row 


1 


0.6956 0.0242 0.6956 






row 


2 


1.3591 0.0362 1.3591 






row 


3 


0.6956 0.0242 0.6956 




Hi 


row 


1 


-1 -1 -2 -1 -1 






row 


2 


-1 -1 -2 -1 -1 






row 


3 


GO 000 






row 


4 


112 11 






row 


5 


112 11 


Diagonal Sobcl 


Ho 


row 


1 


0.3739 0.3022 0.8545 






row 


2 


0.3022 0.6794 0.3022 






row 


3 


0.8545 0.3022 0.3739 




Hi 


row 


1 


-1 -1 -1 -1 0 






row 


2 


-1 -1-10 1 






row 


3 


-1-10 11 






row 


4 


-10 111 






row 


5 


0 1111 
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Once the filter bank is derived, the multiscale MAP 
correlator 20 that works with it can be designed not only to 
detect the feature in the presence of noise, but also to 
distinguish several characteristics of the feature such as 
edges and end-points. Whereas in the past decisions based 
on the outputs of the filter were often heuristic, the multi- 
scale MAP correlator 20 according to the present invention 
is well structuretl and optimized to minimize the detection 
error. 

To reduce the computation burden on the correlator 20 the 
chosen "indicators" in each filtered signal 14a-c are first 
located by the multiscale detector 16. Starting from the 
lowest scale, the detected indicator locations are then traced 
across the different scales. The rule for tracing is to connect 
an indicator location at the lower scale to the spatially 
closest indicator at the next higher scale. This is repeated for 50 
all the observable scales. The result is a "scale-space tree" 
giving us the locations of the indicators at each scale. This 
scale space tree is the input to the MAP correlator 20. 

An optimal and structured MAP correlator 20 is derived 
with Uttle or no information about the signal or the noise. 55 
The decision of the MAP correlator 20 is based on the actual 
filter outputs corresponding to the locations of the indicators 
in the scale-space tree. The correlator 20 is designed to 
minimize the probability of error. To simpUfy the design of 
the MAP correlator, it is assumed that the noise samples are 60 
independent, and identically distributed with a Gaussian 
distribution. The hypothesis testing problem becomes one of 
deciding whether a located feature corresponds to noise or 
an actual signal. Depending on the type of desired 
indicators, the Qorrelator 20 may combine more than one 65 
response from each scale. For example, a zero-crossing is 
marked by checking a sign change in two consecutive 



10 



locations. The two hypotheses are shown to have multivari- 
ate Gaussian distributions. The difference between the two 
densities lies in the mean. The test does not depend on signal 
to noise ratio, but only the desired signal. Details on the 
design of a ID multiscale MAP correlator and a 2D multi- 
scale MAP correlator wiU now be given. 

Assume that the input signal to the multiscale filter bank 
10 has the desired feature embedded in noise. Let us denote 
the input signal samples by x(k), where k represents spatial 
or temporal position. Denote the noise samples an n(k), and 
the feature samples by s(k), 

Jc(A:)=5Cjt)+«(*). 

The noise samples are assumed to be "independent and 
identically distributed" ("i.i.d.") with marginal distribution 
N(0,o^). Therefore, the input samples x(k) are i.i.d. whose 
distribution P can be represented as: 



20 



N{s{k% (T^ ) if signal present 



cr^) i f on] y noi sc present 



It is an i.i.d. sequence with Gaussian distribution. 

The outputs at different scales are Hnear combinations of 
the Gaussian input. Denote the outptit of a filter h. at scale 
i by the samples yj(k). Collect these samples from different 
scales into a vector y;t=[yi(k) y^Qa) . . . Let H represent 
the convolution matrix associated with the filters at different 
scales h„ and the set of samples convolved with H at time 
(position) k, 

H represents a linear transformation appfied to the input. 

The number of input samples needed has to match the 
largest filter. Let n be the size of the largest h, filter, and p 
the number of scales used in the filter. It can be shown that 
the distribution of the linear combination of the n-tuples X;^ 
is a Gaussian random variable. In addition, samples of the 
vector are jointly Gaussian. Their joint density function f 
is characterized by the first two moments. The mean of the 
observation vector y^t can be derived by: 



where m^^ y and m^^, are the means of the input and the 
output respectively. The output correlation matrix R^^ can be 
written in terms of the input correlation matrix R^^^ and H: 

Since the noise samples at the input are i.i.d. with zero mean, 
the covariance matrix is the same as the correlation matrix. 
The mean and the correlation give the density at the output: 



C7) 



( Hst if the £ 
"^*"\0 if only 



feature is present 
noise present 



Denote i(y^ by foCyJ when only noise is present, and denote 
it by fi(yj) in the presence of the desired signal. 

If we denote by Hq the hypothesis when only noise is 
present, and by Hj in the presence of the desired feature, the 
two hypotheses are associated with the densities fg and f^ 
respectively. To determine that a feature is present at time 
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position k, the likelihood ratio test is applied to the two 
densities. L(yjfc)"fi(yfc)/fo(y*)- Assuming uniform costs in 
erroneous detections, the Bayes detector can be reformu- 
lated as: 



Jl if/i 
"\0 othei 



otherwise 



Substituting by the density given above, the MAP correlator 
20 becomes: 



^1 i{{ylR'/^Hs,)i:{l/2)(slH''R'^Hs,) 



(8) 



1^ 0 otherwise 



Note that the test performed by the MAP conelator 20 does 
not depend on the signal to noise ratio, but only on the 
signal/feature sought. This test is useful, when the desired 
feature is indicated by a single position, such as a signal 
maxima. 

The test performed by the MAP correlator 20 is different 
when the featurfe indication is a zero crossing. In this case 
two consecutive positions are indicated for each feature, and 
these must be examined. The two consecutive positions are 
at k and k+1. to give yjt and yj^+i- Tlie two vectors are 
combined by interleaving to give us Yjt«[yi(k)yi(k+l)y2(k) 
y2(k)y2(k+l) - . The indexing represents the different 
scales. Again here, Y^^ can be written as a linear combination 
of the input samples x,^ Recall that y^.=HXj^. Assume that the 
largest filter h, has size n. Take n+j samples of the input 
centered at k, i.e. one extra sample added to Xj^, denote it by 
Let 



h'- 0 

0 H^' 

h^' 0 

0 H^' 



Then 



45 



The MAP correlator 20 for zero crossings can then be 
derived using the above equation (8), where H is replaced by 
H', yk by and x^ by X^. 

Detection of features in 2D signals, such as an image, is 
based on the indications collected across all scales of the 
multiscalc filter 10 and detector 16. In this case, locating the 
indicators also depends on the direction of the feature. For 
example, maxima occur along the gradient direction of 
edges. Once these locations are detected at all the scales, the 
observations at those locations are linked to form the scale 
space tree. This tree is then fed to the MAP correlator 20 to 
decide whether the feature is present or not. 

Assume that the input samples to the filter bank in FIG. 
1 has the desired 2D feanire (signal) embedded in noise. 
Denote those input samples by x(k,l), where k and 1 repre- 
sent the spatial coordinates, the noise samples by n(k,l), and 
the feature samples corresponding to hypothesis j by Sy(k,l), 
we have 

x(AiO-s>(A,0+"CA.O 

The noise samples are assumed to be independent and 
identically distributed ("i.i.d.") with marginal distribution 
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N(0,a^). Consequently, the input samples x(k4) are also Li.d. 
and the distribution P^corresponding to hypothesis j can be 
represented as: 



N (5j (/:,/), cr^) if signjil corresponding to 
hypothesis j preiejit 

/V(0, cr^) if only noise present 



^0 It is an i.i.d. sequence with Gaussian distribution, the 
difference between the hypotheses lies in the mean. 

The filter outputs at different scales are linear combina- 
tions of the Gaussian input. Denote the output of the filter H,- 
at scale I by the samples y^'\k,l). Collecting these samples 

15 from different scales into a vector 

y(k,l)=[y^^^(k,l)y<^>(k,l) . . . Let represent the convo- 
lution matrix associated with the filters H^- from different 
scales, then 



20 



35 



since *H represents a linear transformation applied to the 
input. 

The input samples size matches the size of the largest 
filter Let N be the size of the largest H^- filter, and r be the 
number of scales used in collecting the observations. It can 
be shown that the distribution of the linear combination of 
the n-tuples x(k,l) is a Gaussian random variable, and the 
samples of the vector y(k,l) are jointly Gaussian. Their joint 
density function f is characterized by the first two moments. 
Hie mean of the observation vector y(k,l) is 



where ^y(k,r)i and va^n^^ are the means of the input and the 
output respectively. The output correlation matrix R^^ can be 

written in terms of the input correlation matrix R^. and 7{ : 



Since the noise samples at the input are Li.d. with zero mean, 
the covariance matrix is the same as the correlation matrix. 
The mean and the correlation give the density function of 

yW). 



50 
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im. 0)=- 



(9) 



7/s;(A, /) if signal corre^onding to 
hypothesis y is present 

0 if only noie is present 



Denote f(y(k,l)) by fo(y(k,l)) when only noise is present, and 
denote it by f/y(kl)) in the presence of the signal from 
60 hypothesis j. 

Assuming uniform costs in enoneous detections, the 
Bayes detector can be reformulated as: 



65 



\0 ot 



otherwise 
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With the density given above, the 2D MAP correlator 20 is 



(10) 



0 otherwise 



The hypothesis testing problem becomes one of deciding 
whether a located feature corresponds to noise or actual 
signal. The two hypothesis have multivariate Gaussian dis- 
tributions with different mean. The test does not depend on 
signal to noise ratio, but only on the desired signal. The 2D 
MAP correlator- 20 decides that a feature is present at 
position (k,l) if 



(11) 



15 



20 



25 



yik, if - H-^Hsiik, 0 i -5j {k, l/H'^R-^'Hsafc, 0 



where: 

y(k,l)^ is an indicator location vector at kj; 
Ry^"^HSj(k,l) is a correlation vector at k,l; and 



-isiikJ)fii''R;'y7isiik.O 



is a threshold value. 30 
Note that the test does not depend on the signal to noise 
ratio, but only on the feature signal sought. 

DESCRIPTION OF THE PREFERRED 
EMBODIMENT 



35 



40 



The preferred^ embodiment of the invention is employed 
to detect edges in images produced by an idtrasonic imaging 
system. Referring particularly to FIG. 6 the image to be 
processed is received at input 100 as a 400 by 400 pixel array 
of image intensity data. This image data is applied to one 
input of a convolver 102 that forms the first stage of a 
midtiscale filter bank indicated by dashed line 104. The 
multiscale filter 104 includes three additional convolvers 45 
106, 108 and 110 which are connected in series to form the 
remaining three stages of the multiscale filter 104, Decima- 
tors 112, 114 and 116 connect the successive filter stages and 
sample every other pixel value to change the scale of 
successive filter'stages. The outputs of the convolvers 102, 50 
106, 108 and 110 form the respective multiscale outputs 
120-123, where: 

Scale 1 is at 120; 

Scale 2 is at 121; , 
Scale 3 is at 122; 
Scale 4 is at 123; 

The convolvers 102, 106, 108 and 110 are identical and 
are constructed using digital signal processing ("DSP") 
components. The convolvers perform a convolution opera- 60 
tion on the image array received at one of its inputs with a 
convolution kernel applied to its other input. The convolver 
102 receives a convolution kernel HI from an EPROM 126 
and the remaining convolvers 106, 108 and 110 receive a 
convolution kernel Hq. In this preferred embodiment a Haar 65 
filter function kernel is employed and Hq and Hj are derived 
as described above to yield the following: 
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-1 -1 

-1 -1 

-1 -1 

-1 -1 



It can be appreciated that other filter functions can be 
implemented and their kernels stored in the EPROM 126 for 
application to the multiscale filter 104 when the need arises. 

Each convolver 102, 106, 108, 110 can be implemented 
by Fourier transforms of the image array and the convolu- 
tion kernel applied to its inputs, matrix multiplies them, and 
inverse Fourier transforms the result to produce the filtered 
output image array. 

The deciraators 112, 114 and 116 are also constructed 
using DSP components. They implement the following 
decimation matrix: 



to sample every other pixel as illustrated by the sampling 
lattice in FIG. 5. 

The filtered image arrays produced at the filter outputs 
120-123 are apphed to the inputs of a multiscale detector 
indicated generally at 130. Each of the filtered images is 
processed by a separate detector channel comprised of 
respective detectors 131-134 and masks 136-139. Each 
detector 131-134 scans the filtered image array applied to its 
input and indicates the location therein of the "indicator 
response" of interest. In the preferred embodiment a peak 
detector is employed, and the detectors 131-134 indicate the 
location of peaks, or maxima in the filtered image intensity 
data. These locations are output to the corresponding mask 
136-139, which forms an "indicator response table" by 
reading the filtered image intensity value at each indicated 
"peak" location. Each channel of the multiscale detector 130 
thus produces an indicator response table such as that shown 
at 140 in FIG. 7 which lists the location of each peak and the 
value of the peak. In the preferred embodiment, the func- 
tions performed by the multiscale detector 130 are carried 
out by a general purpose processor under the direction of a 
stored program. 

The indicator response tables 140 are input to a feature 
signature generator 150. In the preferred embodiment the 
function performed by this generator 150 is carried out by a 
stored program which is executed to perform the steps 
indicated by the flow chart in FIG. 9. This program examines 
the indicator response table 140 for each scale, and produces 
a scale space tree ("SST") table which is illustrated at 142 
in FIG. 8. 

Referring to FIG. 9, the first step in the procedure is to 
initialize the relevant data structures as indicated by process 
block 160. This includes, for example, clearing the SST table 
142. The indicator responses 140 for the first scale are then 
stored in the first row of the SST table 142 as indicated at 
process block 162. This is the highest resolution scale and 
the indicator response table 140 for scale 1 will contain the 
most entries. A scale pointer is then set to read the indicator 
response table 140 for scale 2, as indicated at process block 
164, and a loop is entered in which the indicated peaks in 
scales 2, 3 and 4 are examined. 
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Referring still to FIG. 9 each indicator response table 140 
is examined by setting a response pointer to "1" at process 
block 166 to read the location of the first indicated peak. If 
this location corresponds (plus or minus one pixel location) 
with a peak location indicated in the previous scale, as s 
determined at decision block 168, the value for that location 
and scale is loaded into the SST table 142 as indicated at 
process block 170. Otherwise, the system loops back to 
examine the next indicator response after incrementing the 
response pointer at 172. After the last indicator response in lo 
the table 140 has been examined, as determined at decision 
block 174, the scale pointer is incremented at 176 and the 
indicator response table 140 for the next scale is examined. 
The process terminates when the last scale is finished as 
determined at decision block 178. 15 

Referring particularly to FIG. 8, the resulting SST table 
142 stores a set of n pixel locations at which peaks were 
found in the scale I indicator response table 140. Where 
corresponding peaks were found at higher scales, values are 
entered for those scales. It should be apparent that locations 
where peaks were found at all scales have a higher prob- 
ability of being the location of an edge, whereas locations 
with indicator values at only scale I arc likely noise or 
weaker edges. The values at each SST table location (i.e. 
each column in FIG. 8) is referred to herein as an "indicator 
location vector", and it is this vector which is examined to 
determine if an edge is present at that location. The values 
in each indicator location vector can be viewed as a "sig- 
nature" of some feature in the image produced by the 
multiscale filter J04. 

Referring again to FIG. 6 the SST table 142 produced by 
the feature signature generator 150 is examined by a MAP 
correlator 180 to determine which indicator location vectors 
therein represent true edges. The functions performed by the 
MAP correlator 180 are carried out by a general purpose 
computer which executes a program illustrated by the flow 
chart in FIG. 10. 

Referring particularly to FIG. 10, the first step performed 
by the MAP correlator 180 is to initialize the data stmctures 
as indicated by process block 182. Such structures include, 
for example, a feature image array which is set to zero, or 
"black". As indicated by the above equation 11, a correlation 
vector is then calculated at block 184. This correlation 
vector may be viewed as the ideal indicator location vector 
produced by the multiscale filter 104 in response to an edge 
in the input image. Next, a threshold value is calculated at 
process block 186. This calculation is also indicated in the 
above equation 11. A loop is then entered in which each 
indicator location vector in the SST table 142 is read at 
process block 188, cross correlated with the correlation 50 
vector at process block 190, and the calculated correlation 
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value compared with the threshold at decision block 192. If 
the indicator location vector is similar enough to the corre- 
lation vector that the calculated correlation value exceeds 
the threshold, an edge is present and the pixel at the vector 
location is set to a bright value at process block 194. 
Otherwise, process block 194 is skipped. This process 
continues until each indicator location vector in the SST 
table 142 is examined. When the last vector has been 
examined as determined at decision block 196, the feature 
image is complete and it is output at process block 198 to a 
display 200 (FIG. 6). Hie feature image indicates each 
detected edge as a bright line. This may be the final 
objective, or further steps may be carried out using this 
feature image to locate and identify shapes, and calculate 
their areas or volumes. 
What is claimed is: 

1. A detector system for locating a feature in an input 
signal, the combination comprising: 

a multiscale filter which applies a filter operator to the 
input signal to produce a plurality of filtered output 
signals at a corresponding plurality of scales; 

a multiscale detector which receives the filtered output 
signals and produces a table which indicates the pos- 
sible locations in the input signal of the delected feature 
therein; and 

a MAP correlator which examines the table and correlates 
the values therein at each indicated feature location 
with a model response, said MAP correlator producing 
an output signal indicating the location of the feature if 
the correlation exceeds a preset threshold. 

2. The detector system as recited in claim 1 in which the 
input signal is a two-dimensional image and the output 
signals form a bit map. 

3. The detector system as recited in claim 1 in which the 
multiscale detector includes a feature signature generator 
which produces the table by forming a scale space tree from 
the corresponding locations in each filtered output signal at 
which the feature is detected. 

4. The detector system as recited in claim 1 in which the 
multiscale filter is a filter bank comprised of a set of series 
connected filter stages. 

5. The detector system as recited in claim 4 in which each 
filter stage includes a convolver which receives the input 
signal and convolves it with a filter kernel, and in which the 
first filter bank stage employs a filter kernel and the 
remaining filter bank stages employ a filter kernel Hq. 

6. The detector system as recited in claim 1 in which the 
input signal is a series of two-dimensional images and the 
output signals form a scries of bit maps. 
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